Shock statistics in higher- dimensional Burgers turbulence 
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We conjecture the exact shock statistics in the inviscid decaying Burgers equation in D > 1 
dimensions, with a special class of correlated initial velocities, which reduce to Brownian for D — 
1. The prediction is based on a field-theory argument, and receives support from our numerical 
calculations. We find that, along any given direction, shocks sizes and locations are uncorrelated. 

PACS numbers: 



Decaying turbulence is characterized by the existence 
of an inertial range in the inviscid limit (small viscosity 
limit) with scaling and multi-scaling. These features are 
shared by simpler models as passive advection [1] and the 
Burgers equation [2] where the velocity field present man- 
ifolds of discontinuities called shocks. How the velocity 
field evolves from a prescribed random initial condition 
v(r,t = 0), f € MP , and the statistics of these shocks, 
are of high interest. Unfortunately, exact results are re- 
stricted to few solvable cases in space dimension D = 1 
[2"HTT] or in the limit D = oo [TH |T3] . In this Letter we 
present a solution for generic D, for a non-trivial class 
of random initial conditions. These appear naturally in 
related works in the context of elastic interfaces in D = I 
[MlflT] and more recently in higher D [17] . At this stage 
the solution is a conjecture, based on a field-theory ar- 
gument. Here, we state the conjecture and provide an 
accurate numerical test. 

The decaying Burgers equation describing a potential 
flow-velocity field v(r, t) — VV(r, t) reads 



d t v = vV 2 v W 2 

2 



(1) 



The inviscid limit corresponds to v = + [18] . Our model 
is defined by choosing the distribution of the initial ve- 
locity field v(f, 0) as a centered gaussian with stationary 
increments 5v(fi,?2) = v(fi, 0) — v(r2,0) of correlations: 
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B, 



6vi(f ,f +r)Svj(fo,f + r) = -^\r\(5ij+fifj) (2) 

with statistical transnational invariance, where f = r/H- 
We denote by 7TT averages over initial conditions. In 
D = 1 this reduces to a Brownian initial velocity [3HS]- 
In general D there is no obvious Markov property, except 
that the velocity along any given direction is a Brownian. 

We now state our prediction, obtained from field- 
theoretical considerations, and to be considered a conjec- 
ture in the mathematical sense: The characteristic func- 
tion of the 2-point distribution of velocity increments is 



-X-[iT(a:Si,t)-ff(0,t)] = p^^fAJ-AJ 



(3) 



In the notations of equation ^ we set without loss of 
generality r = 0, and r = xe\, x > 0; e\ is the unit 



vector along the x-axes. The generating function Z t (X) 
is a function of two arguments, Z t (X) — Z t (X x , X]_), 
X± = X — X x ei, and will be discussed below. Further- 
more, Z t (X) = S m Z(S m X/t) where S m = Bt 2 is a spatial 
scale related to a characteristic size of the shocks defined 
below. The function Z(X) is given below in some spe- 
cial directions, its full expression, not reproduced here, 
is computed in |17j . The general form ^ can be derived 
under the following assumptions : (i) velocity increments 
are localized in shocks; (ii) along a given direction, the lo- 
cation of shocks are independent (i.e. Poisson statistics) 
(iii) shock sizes are mutually uncorrelated and indepen- 
dent from locations. The converse is also true [T5] . 

Both S rn and Z(X) are related to shocks as follows. Let 
us denote by x a the discrete set of points where the shock 
manifold intersects the x axis, and (minus) their associ- 
ated velocity jumps as Si = t[v(x a — + ) — v(x a + + )] 
where a labels the shocks. The shock-size density 

is p(S) = J2 a 3( x ~ x a)S(S — S a ). The characteristic 
shock size S m is defined from the moments of the longi- 
tudinal shock component S x as 

(S 2 X ) 



2(S X 



(4) 



where (• • • ) denotes averages over the shock density p(S) 
which takes the form 



P0) = -^-p(q- 



(5) 



p(s) is a function of the reduced shock-size s :— S/S m . 
By construction (s 2 ) = 2. The identity / dS p(S)S x = 1 
implies a second normalization condition (s x ) = 1. Note 
that here and below (. . .) denote moments either over 
p(S) or p(s). Expanding Eq. ([3]) for small x and using 
(|5| one finds that Z(X) is the generating function for the 
distribution of reduced shock sizes, 

Z{X) := - l) := / ds ? - l) p(s) . (6) 

Its expansion has been computed in |17j : 

Z(X) = A K +iA 2 + iA 2 + 2A K A 2 +^(A 2 ) 2 +^ 2 A 2 -A^+. . . 

(7) 
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It implies universal moment ratios, in particular 

(Si) 2 2 



(sr. 



D - 1 



(8) 



where here S± and s± denote the component of the shock 
orthogonal to x. While the set of shocks along x are mi- 
correlated both in position and size, a property which 
indeed implies by contrast, longitudinal and trans- 
verse components of a given shock are correlated, as from 
([7]) one can calculate higher moments, e.g. 



(S x Sj_) (S x 
(SI) 2 



= 4(D-1) . 



(9) 



We now indicate the origin of our conjecture, by recalling 
the connection to disordered systems. Eq. ([!]) is solved 
by the Cole-Hopf transformation [2] in the limit v —> 0: 



V(f,t) = 



2t 



(u - rf + V(u) 



(10) 



where V(u) is the potential associated with the initial 
condition, i.e. v(r,t = 0) = W(r). Hence for a ran- 
dom initial condition the problem is equivalent to find- 
ing the minimum energy position of a particle in a ran- 
dom potential, plus a harmonic well. Denoting by u(r) 



the position of the minimum in (10), the velocity field is 
v(f,t) = [f—u(r)]/t. At the shocks, the minimum jumps, 
and the shock size is S — u(x a + + ) — u(x a — + ). Note 
that u(r) = r which implies (S x ) — 1 as stated above. 

The random potential V(u) corresponding to the 
present model ^ is a generalization of the ID random 
acceleration process PHI HI] to D dimensions. To define it 
one needs a large-scale regularization; we choose periodic 
boundary conditions of period L in all D directions, 



V{u) = IT 



qjtO 



'-+H 



(11) 



where q = f-ft, n e {-£/2 + 1, . . . , i/2 - 1, i/2} , 
in the limit L — > 00, and H = 3/2. In real space 
this leads to a non-analytic cubic potential correlator 
V{u)V{u') = R (u-u") with R Q (u)~R {0) 



\A L u 2 



f M 3 + 0(1/L) with A L = 0.0182Lcr 2 + O(L ) and 
B = <7 2 /(37r) + 0(1/ L). The initial velocity correlator is 
Vi(r,t = 0)vj(0, t = 0) = —didjRo(r) with independent 
increments distributed as in ([2]). 

In a nutshell the basis for the conjecture is as follows. 
The present model is the d — limit of a model of an 
elastic manifold (of internal dimension d) in a random 
potential and a quadratic well of curvature 1 jt. The anal- 
ogous variable to u(f) is the center of mass of the mani- 
fold, and V{r) the energy of the optimal configuration as 
a function of well position. Its second cumulant defines 
a renormalized potential disorder correlator R(r) for any 
d, which is shown to obey a Functional RG equation as 




FIG. 1: Left. Step One: reduction to one dimension. Effective 
1-dimensional potential V(i, i) after the minimization over j, 
as given in Eq.( |14[ ), at different times (from top to bottom 
t = 0,2, 8). The potential becomes deeper and deeper as time 
increases. Right. Step Two: location of the minimum. Solid 
stair-case line is i m in(x), dashed line is jmin(x) for t = 8. The 
drift x is indicated. Shocks are only forward in x direction. 



t is varied. This equation can be solved perturbatively 
ini?inad = 4 — e expansion. It turns out that the 
initial correlator Ro(r) corresponding to ^ solves the 
FRG equation to all orders in e, i.e. there are no loop 
corrections. This implies that the correlation functions 
need only be computed to tree-level, either by recursion 
or from a saddle-point method, as detailed in |17j . This 
leads to ^ and to Z t (X), which hold for any D and any 
d, for this choice of initial conditions, although we need 
only d = (Burgers). A further result, proved to low- 
est order in e = 4 — d [T7] but which we expect to hold 
for any d, is that ([2]) is an attractive fixed point of the 
RG, hence for velocity correlations which differ from ^ 
only at small r, the behaviour at large t again follows 
([3]) [IS] . Of course we cannot exclude non-perturbative 
corrections, hence our prediction is, strictly, a conjecture. 
In support we note that for D = 1 it has been proven in 
[6]. To check it in D = 2 we now turn to numerics. 

A powerful algorithm allows to solve this problem for 



a slightly modified version of Eq. ( 10 ), with a discretized 
variable u — and a continuous variable r — xe\ 



V(xei,t) = min 

l<i,j<L 



(i - xf 

2~r~ 



2t 



V(i,j) 



, (12) 



for any x in the interval (0, L). Let us now discuss how 
the algorithm finds the site ■u m i n (a;) = {imin{%) > jmin(%)) 
which satisfies the minimization condition (12): 



Step 1: Reduction to a 1-dimensional problem. For 
each value of i we perform a minimization over the trans- 
verse coordinate j, keeping in memory the location of the 
minimum, j^ in (i). Since this operation does not involve 
x, the effective dimension of the problem is reduced to 1, 
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FIG. 2: Color Online. Left : Convergence of the measured 
Z™(J) to the analytical prediction (16 1 for different times: 
t = 2.04 (filled circles), t = 2.78 (open circles), t = 4 
(squares). Right: measured Z 2{X) for t = 4 (squares) com- 
pared to the prediction \17\ (solid line). 



and Eq. ( 12 ) becomes 



V(xei,t) = min 

Ki<L 



V(i,t) = min 

l<j<L 



{i~xf 



2t 



+ V(i,t) 



2t 



V(i,j) 



(13) 
(14) 



The reduced potential V(i, t) is plotted on Fig. [T] (left). 

Step 2: Determination of i m i n (x). The latter is an 
increasing piecewise constant function of x. The mini- 
mum location in the original D = 2 lattice is given by 
jminO) = j,nin( i min(a;)) : For x = the minimum posi- 



tion is found from Eq. (13 1. Increasing x, the minimum 
remains in imi n (x = 0) up to a threshold xi, above which 
the minimum takes a new value i m i n (xi) > «min(0). For 
all i > i m j n (x = 0) we find the value of x satisfying 



2t 



■V(i,t) 



(»min(0) - z) 

2i 



+ V(w(0),t) . (15) 



xi is the smallest value of x for which this condition is 
satisfied. One then searches the next minimum and the 
procedure is iterated up to x — L, see Fig. 1 (right). 

Step 3: Shock sizes. Given the sequence of minima 
locations u m in(x) = (i m in(x), j m in{x)), the shocks sizes 
S are the discontinuities in these piecewise functions of 
x. The velocity profile is v x (x,t) = [x — i m in(x))/t, 
v y (x,t) = —jmm( x )/t- Note that in our discrete model 
the shock size is cut off from below at So = 1 an d from 
above at L. Self-affine scaling and Eq. ([I]) are expected 
to hold in the continuum limit when So <C S m < L or 
equivalents S /S m = l/{Bt 2 ) < s < L/(Bt 2 ). 

(iv) Step 4: Numerical implementation. In practice, we 
consider a D = 2 square latice (usually of size L = 2 12 ), 
the correlated random potential V(i,j) is constructed 
from L 2 independently distributed Gaussian random 



numbers via a "fast Fourier transform" of Eq. (111. Note 




FIG. 3: Color Online. Left. Plot of s 2 x pi(s x ) (top curve ) 
and s 2 L [p2(s±) +P2(— s±)] (bottom curve). Shocks at t = 8 
(S m = 7 ± 0.1). Solid line represents the analytical predic- 
tions. Right. Test of the ratio Squares are for L = 2 8 , 
circles for L = 2 10 , and triangles for L = 2 12 . Numeri- 
cal data approach the predicted value when L is large. For 
L/Sm — > the particle feels the periodic potential, i.e. S x = L 
and S± = 0. For t — > 0, lattice spacing is important: S x = 1 
and Sj_ = 0. The plateau value is consistent with 4 — cil/ _1//2 . 



over integers from — L/2 + 1 to L/2. The zero mode n = 
is set to zero, Vq = 0. We choose a 2 = 1, which implies 
that B = l/(37r) in formula pj. We collected a large 
number of shocks (~ 10 6 — 1(F) using many samples, 
from which we computed S m and verified the prediction 
S m = t 2 /(3ir). From the reduced sizes s a := S a /S m , 
we measured Z(X X ,X±) = jjJ2 a ( eX Sa ~~ 1)> specifically 
Zi(A) := Z(\,0) and Z 2 (X) := Z(0,X). The conjecture 
states that for the longitudinal component of the shock 
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1 



(16) 



with pi (s x ) :— J ds±p(s Xl s±), i.e. the same as obtained 
for D = 1 H [TU] and for the related Galton process 
[161 122]- This is verified on Fig. [| (left) and Fig. [| (left). 
Since the agreement is very good, we have plotted on 
Fig. [2] (left) the difference with the analytical prediction 
to emphasize the small deviations. These deviations are 
more important for large negative A ~ — 1/so, sensitive 
to the small lattice cutoff so = So/S m for the reduced 
shock sizes. Increasing the time, sq decreases as sq ~ 
1/t 2 . The prediction for the characteristic function of 
the y component of the shock sizes, Z 2 (A), is obtained 
by elimination of 9 in the system of equations 



X(6) 



Z 2 (0) 



sin 6- 



y/5 - cos(4<9) + 2 



[1 - cos(2(9) + ^5 - cos(46>)] ' 
OS 9 ^5 - cos(46>) - 2 
2 1 - cos(2(9) + y/5 - cos(4<9) ' 



(17) 



(18) 



that the sum over the components of ft are now running 



Numerically, the Laplace inversion can be performed 
to determine p 2 (s±) = j ds x p(s x , s±) with high pre- 
cision. (It is an integral over a segment of 9 in the 
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complex plane.) P2(s) is plotted on Fig. [3] (left). For 
large s, p 2 {s) « 1.7304|s|- 5 / 2 e-°- 2698 l s l ; while for small s 
P2{\s\) = 0.12375|s|~ 3 / 2 . We have plotted the measured 
and calculated Z 2 {\) on Fig. || (right). Since p 2 (s) is sym- 
metric in s, the same holds true for Z 2 (X). The left and 
right edges of the analytic curve are at |A*| = 0.2698..., 
the constant in the exponential decay of pi{s). The agree- 
ment is excellent up to this point, where the size L cuts 
the divergence for |A| > A*. 

We now discuss shock correlations: First, the universal 
ratio ^ was measured to be 2.034 ± 0.015, very close to 
its analytical prediction. Second, correlations of jumps 
in the different directions are measured by the ratio 
plotted on Fig. [3] (right). In both cases, the deviations 
can be attributed to finite-size corrections, see the cap- 
tion of Fig. 3. Third, we studied the correlations between 
subsequent shocks. To emphasize the remarkable nature 
of the present model (H = 3/2) we compare with two 
other ones, with potential given by (11) with H = 0.5 
and H = 1. In Fig. (top) we show the connected cor- 
relation (s Xia s Xia + p ) c of the longitudinal size s x>a of a 
shock with the p-th subsequent shock. Fig. [1] (bottom) 
focuses on the correlations between the location and the 
size of shocks. For example, large shocks are more iso- 
lated with respect to small shocks? To check this we 
compute the average distance {x a +i — x a ) between con- 
secutive shocks, normalized by its averaged and called 
A(s x ), as a function of the longitudinal size s x of the 
shock a. Fig. [4] (bottom) shows strong correlations for 
H = 0.5 and H = 1. For H = 3/2 no effect was detected. 
This remarkable statistical independence of the shocks is 
essential for the main formula 

To conclude, we proposed a conjecture for the D- 
dimensional decaying Burgers equation with initial con- 
ditions which generalize the Brownian for D = 1. We 
tested it numerically and checked that the shocks and 
the velocity increments along an axis are statistically in- 
dependent at any time t. The conjecture is based on van- 
ishing loop corrections in the field theory for the disor- 
dered problem and its generalization to elastic manifolds. 
Although confirmed within our numerical accuracy, any 
deviation would have important consequences for the - 
probably non-perturbative - corrections to the field the- 
ory. We hope this motivates efforts to prove or infirm 
our conjecture on a rigorous basis. 
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